Compare datasets

---

library(raster)

library(plyr)
library(dplyr)
library(reshape)

library(ggrepel)
library(ggplot2)
library(ggpmisc)
library(ggpubr)

library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "NA" = "#D5D8DC")

col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

1 Data

Liens utiles:

• https://www.fao.org/3/x0490E/x0490e0k.htm#annex%203.%20background%20on%20physical%20parameters%20used%20in%20evapotranspiration%20computatio

# era5vA <- read.table("era5_AI.txt") %>%
#   mutate(tas = tas-273.15, pr = mpr * 60 * 60 * 24 * 365, rsds = ssrd / (60*60*24) * 1e-6 * 60*60*24) %>% # conversion from J/m2 to J/sm2 to MJ/m2/day  
#   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","rsds","ea","ET0","AI","cat.AI"))
# 
# 
# write.table(era5vA, "era5vA.txt")

era5vA <- read.table("ERA5/era5_AI.txt") %>% 
  mutate(z = z, tas = t2m - 273.15, sfcWind = wind, model = "historical", source = "ERA5") %>%  #pr, ET0 already in mm/year
  select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)

write.table(era5vA, "era5vA.txt")

wdcvA <- read.table("wdc_AI.txt") %>% mutate(z = z, ea = vapr, ET0 = ET0) %>% # srad in MJ/m2/day, pr and ETO already in mm/year
 select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI","cat.AI")) %>%
  filter(lm == 1)


write.table(wdcvA, "wdcvA.txt")


cmip6 <- read.table("cmip6_AI.txt")

cmip6vA <- cmip6 %>% subset(period == "1970_2000") %>%
   group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
   dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(spr, na.rm = T), # pr in mm/year
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>% # ET0 in mm/year
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
   select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

cmip6vA$cat.AI <- cmip6vA$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6vA$cat.AI[which(cmip6vA$ET0 < 400)] <- "Cold"

write.table(cmip6vA, "cmip6vA.txt")

era5vA <- read.table("figures_article/era5vA.txt")
wdcvA <- read.table("figures_article/wdcvA.txt")
cmip6vA <- read.table("figures_article/cmip6vA.txt")

2 Regions and Continents

colors_continents <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
                  "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
                  "#aec7e8", "#ffbb78", "#98df8a", "#ff9896")

ggplot(era5vA)+geom_raster(aes(x= lon, y = lat, fill = Continent))+
  scale_fill_discrete(type = colors_continents)+
  theme_void()

“Continents” to remove are: * Arctic (too little gridcells) * Atlantic (ocean) * Indian (ocean) * Pacific (ocean) * Southern (ocean)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
plot(land_mask)

land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 

ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

ipcc_sf <- st_as_sf(ipcc_regions)

ipcc_sf$Name2 <- ipcc_sf$Name %>% str_wrap(width = 10)

ggplot(data = filter(ipcc_sf, Type %in% c("Land","Land-Ocean")))+
  geom_sf(aes(fill=Continent))+
  scale_fill_manual(values = colorvec)+
 geom_sf_text(aes(label = Name), size = 3, position = position_jitter(height = 3, width = 3, seed = 1))+
  #geom_sf_text(aes(label = Acronym), size = 2)+
  theme_void()+
  theme(legend.position = "bottom")

#ggsave("ipcc_regions.png", width = 12, height = 8, units = "cm", dpi = "retina")

3 Aridity Index FAO

3.1 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with country borders

calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
#col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"

calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "ERA5, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top"),

nrows = 4, ncol = 1

)

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("worldclim_AIref.png", width = 12, height = 8, dpi = "retina")

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("era5_AIref.png", width = 12, height = 8, dpi = "retina")

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    borders(colour= "grey60")+
    labs(title = "", fill = "")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggsave("cmip6_AIref.png", width = 12, height = 8, dpi = "retina")

wdcfordiff <- calib %>% subset(source == "Worldclim") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.wdc")) 
era5fordiff <- calib %>% subset(source == "ERA5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.era5")) 
cmip6fordiff <- calib %>% subset(source == "CMIP6") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmip6")) 

diffcatAI <- wdcfordiff %>% merge(era5fordiff, by = c("lon","lat")) %>% merge(cmip6fordiff, by = c("lon","lat"))

diffcatAI$era5_wdc <- with(diffcatAI, 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.wdc %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.wdc == "Hyperarid", "dryer",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.wdc %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.wdc %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.wdc %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.wdc %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.wdc == "Cold","wetter",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.wdc %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Cold" & cat.AI.wdc %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Hyperarid" & cat.AI.wdc %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.era5 == cat.AI.wdc, "same",
          NA))))))))))))

diffcatAI$era5_cmip6 <- with(diffcatAI, 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.era5 == "Arid" & cat.AI.cmip6 == "Hyperarid", "dryer",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.era5 == "Semi-arid" & cat.AI.cmip6 %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.cmip6 %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.era5 == "Dry subhumid" & cat.AI.cmip6 %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.cmip6 == "Cold","wetter",
          ifelse(cat.AI.era5 == "Humid" & cat.AI.cmip6 %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Cold" & cat.AI.cmip6 %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.era5 == "Hyperarid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.era5 == cat.AI.cmip6, "same",
          NA))))))))))))

diffcatAI$wdc_cmip6 <- with(diffcatAI, 
          ifelse(cat.AI.wdc == "Arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.wdc == "Arid" & cat.AI.cmip6 == "Hyperarid", "dryer",
          ifelse(cat.AI.wdc == "Semi-arid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.wdc == "Semi-arid" & cat.AI.cmip6 %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.wdc == "Dry subhumid" & cat.AI.cmip6 %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.wdc == "Dry subhumid" & cat.AI.cmip6 %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.wdc == "Humid" & cat.AI.cmip6 == "Cold","wetter",
          ifelse(cat.AI.wdc == "Humid" & cat.AI.cmip6 %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.wdc == "Cold" & cat.AI.cmip6 %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.wdc == "Hyperarid" & cat.AI.cmip6 %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.wdc == cat.AI.cmip6, "same",
          NA))))))))))))

ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.wdc)) + geom_raster(aes(x=lon, y = lat,  fill = era5_wdc))+
    scale_fill_manual(values = c("same" = "white", "dryer" = "#D73027", "wetter" = "#DAE2F7"), na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-Worldclim, 1970-2000", fill = "Worldclim dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = era5_cmip6))+
    scale_fill_manual(values = c("same" = "white", "dryer" = "#D73027", "wetter" = "#DAE2F7"), na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.wdc != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = wdc_cmip6))+
    scale_fill_manual(values = c("same" = "white", "dryer" = "#D73027", "wetter" = "#DAE2F7"), na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories Worldclim-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to Worldclim")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.wdc)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.wdc))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-Worldclim, 1970-2000", fill = "Worldclim dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.era5 != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.cmip6))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories ERA5-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to ERA5")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

ggplot(filter(diffcatAI, cat.AI.wdc != cat.AI.cmip6)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI.cmip6))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    borders()+
    labs(title = "Difference aridity categories Worldclim-CMIP6, 1970-2000", fill = "CMIP6 dryer or wetter compared to Worldclim")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "bottom")

3.1.1 Compare ETO and pr in WDC, ERA5, CMIP6

bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.pr <- c("(-400,0]" = colvecrdbu[5], "(0,400]"= colvecrdbu[6],"(400,800]"= colvecrdbu[7], "(800,1.2e+03]" = colvecrdbu[8], "(1.2e+03,1.6e+03]"= colvecrdbu[9], "(1.6e+03,5e+03]" = colvecrdbu[10])

calib$pr.cat <- calib$pr %>% cut(breaks = bpr) 

g.pr.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation in Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
      scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation in ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.pr.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual precipitation CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.pr <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = pr.cat))+
    scale_fill_manual(values= col.pr, na.translate = F)+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))

#)

bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <-  brewer.pal(n = 11, name ="RdYlBu")

col.et0 <- c("(-400,0]" = colvecrdbu[8], "(0,400]"= colvecrdbu[7],"(400,800]"= colvecrdbu[6], "(800,1.2e+03]" = colvecrdbu[5], "(1.2e+03,1.6e+03]"= colvecrdbu[4], "(1.6e+03,5e+03]" = colvecrdbu[3])

calib$et0.cat <- calib$ET0 %>% cut(breaks = bpr) 

g.et0.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration Wordclim", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration ERA5", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

g.et0.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    borders(colour= "grey60")+
    labs(title = "Annual evapotranspiration CMIP6", fill = "mm/day")+
    ylim(-55,90)+
    theme_void()+
    theme(legend.position = "none")

l.et0 <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = et0.cat))+
   scale_fill_manual(values= col.et0, na.translate = F)+
    geom_raster(aes(x = lon, y = lat), fill = "white")+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "top", legend.key.width = unit(2, "cm"))

ggarrange(
  ggarrange(plotlist = list(g.pr.wdc,g.pr.era5, g.pr.cmip6, l.pr), nrow = 4, ncol = 1),
  ggarrange(plotlist = list(g.et0.wdc,g.et0.era5, g.et0.cmip6, l.et0), nrow = 4, ncol = 1),
  ncol = 2)

3.2 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with Antarctica

3.2.1 Map aridity categories for the 3 databases

calib <- rbind(cmip6vA, era5vA, wdcvA)

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
# col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
col.cat <- c("Hyperarid" = "#D73027", "Arid"= "#FDAE61","Semi-arid"= "#FEE08B", "Dry subhumid" = "#ABD9E9", "Humid"= "#4575B4", "Cold" = "#DAE2F7")

calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"
calib$cat.AI <- factor(calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

ggarrange(
  
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "Wordclim, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "ERA5, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat)+
    labs(title = "CMIP6, 1970-2000", fill = "")+
    theme_void()+
    theme(legend.position = "none"),

ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
  geom_raster(aes(x = lon, y = lat), fill = "white")+
    scale_fill_manual(values = col.cat)+
    labs(title = "", fill = "")+
    theme_void()+
    theme(legend.position = "bottom"),

nrows = 4, ncol = 1

)

3.2.2 Proportion aridity category by subregion, continent et database

calib$same.cat <- NA
calib$cat.ref <- NA

for(i in 1:nrow(calib)){
  loni = calib$lon[i]
  lati = calib$lat[i]
  catrefi <- subset(calib, lon == loni & lat == lati & source == "CMIP6")$cat.AI %>% as.character()
  cati <- calib$cat.AI[i] %>% as.character()
  calib$same.cat[i] <- ifelse(cati == catrefi, "y",
                              ifelse(cati!= catrefi, "n", 
                                     NA))
  calib$cat.ref[i] <- catrefi
}

write.table(calib, "figures_article/calib.txt")

# table_calib_c <- calib %>% select(c("Continent","source","cat.AI")) %>%
#   reshape2::melt(id.vars = c("Continent","source")) %>%
#   group_by(Continent,source,value) %>%
#   summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
#   ungroup()
# 
# table_calib_c$value[is.na(table_calib_c$value) == T] <- "NA"
# 
# ggplot(table_calib_c)+
#   geom_col(aes(x = Continent, y = perc, fill = value, position = "fill"))+
#   scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
#   theme_minimal()+
#   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
#   facet_grid(rows = vars(source))+
#   labs(x = "Continent", y = "%", fill = "Aridity\ncategory", y = "Proportion of aridity category")
# 


table_calib <- calib %>% select(c("Continent","Name","Acronym","source","cat.AI")) %>%
  reshape2::melt(id.vars = c("Continent","Name","Acronym","source")) %>%
  filter(!is.na(value)) %>%
  group_by(Continent,Name,Acronym,source,value) %>%
  summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
  ungroup() %>% filter(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC"))

table_calib$value[is.na(table_calib$value) == T] <- "NA"
table_calib$Continent[which(table_calib$Continent == "EUROPE-AFRICA")] <- "EUROPE"
table_calib$Continent[which(table_calib$Continent == "CENTRAL-AMERICA")] <- "CENTRAL\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "NORTH-AMERICA")] <- "NORTH\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "SOUTH-AMERICA")] <- "SOUTH\nAMERICA"

table_calib$Continent <- factor(table_calib$Continent, levels = c("AFRICA","EUROPE","ASIA","OCEANIA","NORTH\nAMERICA","CENTRAL\nAMERICA","SOUTH\nAMERICA","POLAR","NA"))

table_calib$value <- factor(table_calib$value, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = perc, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "%", fill = "Aridity\ncategory")

ggplot(table_calib)+
  geom_col(aes(x = Acronym, y = count, fill = value, position = "fill"))+
  scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))+
  scale_y_continuous(position = "right")+
  theme_minimal(base_size = 12)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5), 
        legend.position = "bottom", 
        strip.text.x = element_text(size = 6))+
  facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
  labs(x = "", y = "Count", fill = "Aridity\ncategory")

3.2.3 Pie chart

library(ggrepel)

pie_calib <- table_calib %>% group_by(source, value) %>% summarise(gridcells = sum(count)) %>% mutate(percent = round(gridcells/sum(gridcells)*100, 1), lab.y = (rev(cumsum(rev(percent)))) - percent*0.5)


ggplot(pie_calib)+geom_bar(aes(x = "", y = percent, fill = value), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = percent), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

pie_calib <- rbind(cmip6vA, era5vA, wdcvA)

pie_calib$cat.AI <- factor(pie_calib$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold")) 

pie_calib <- pie_calib %>% 
  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("source","cat.AI")) %>%
   table() %>% as.data.frame() %>%
   group_by(source) %>%
    mutate(perc = round(Freq/sum(Freq)*100, 1),
         lab.y = rev(cumsum(rev(Freq))) - Freq*0.5) 

ggplot(pie_calib)+geom_bar(aes(x = "", y = Freq, fill = cat.AI), width = 1, stat = "identity")+
  coord_polar("y", start = 0)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  facet_grid(cols = vars(source))+
  geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
  labs(fill = "")+
  theme_void()+theme(legend.position = "bottom")

# cmip6 <- read.table("cmip6_AI.txt")
# 
# cmip6.ref <- cmip6 %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","source","AI","spr","t","sET0")) %>% setNames(c("lon","lat","model","period","source","AI.ref","spr.ref","t.ref","ET0.ref"))
# cmip6.ref$cat.AI <- cmip6.ref$AI.ref %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
# cmip6.ref$cat.AI[which(cmip6.ref$ET0.ref < 400)] <- "Cold"
# 
# write.table(cmip6.ref, "figures_article/cmip6.ref.txt")

cmip6.ref <- read.table("figures_article/cmip6.ref.txt")

wdc <- filter(calib, source == "Worldclim") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.wdc"))
era5 <- filter(calib, source == "ERA5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.era5"))
cmip6.cal <- filter(calib, source == "CMIP6") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmip6"))

casesm2 <- filter(cmip6.ref, source == "CAS-ESM2") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.casesm2"))
cesm <- filter(cmip6.ref, source == "CESM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cesm"))
cmcc <- filter(cmip6.ref, source == "CMCC") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmcc"))
cmccesm2 <- filter(cmip6.ref, source == "CMCC-ESM2") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cmccesm2"))
cnrm <- filter(cmip6.ref, source == "CNRM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.cnrm"))
ecearth3 <- filter(cmip6.ref, source == "EC-Earth3") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.ecearth3"))
fgoals <- filter(cmip6.ref, source == "FGOALS") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.fgoals"))
gfdlesm4 <- filter(cmip6.ref, source == "GFDL-ESM4") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.gfdlesm4"))
inm <- filter(cmip6.ref, source == "INM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.inm"))
inmcm5 <- filter(cmip6.ref, source == "INM-CM5") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.inmcm5"))
mpi <- filter(cmip6.ref, source == "MPI") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.mpi"))
mri <- filter(cmip6.ref, source == "MRI") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.mri"))
noresm <- filter(cmip6.ref, source == "NorESM-2-MM") %>% select(c("lon","lat","cat.AI")) %>% setNames(c("lon","lat","cat.AI.noresm"))

all <- merge(wdc, era5, by = c("lon","lat")) %>% 
  merge(cmip6.cal, by = c("lon","lat")) %>% 
  merge(casesm2, by = c("lon","lat")) %>%
  merge(cesm, by = c("lon","lat")) %>%
  merge(cmcc, by = c("lon","lat")) %>%
  merge(cmccesm2, by = c("lon","lat")) %>%
  merge(cnrm, by = c("lon","lat")) %>%
  merge(ecearth3, by = c("lon","lat")) %>%
  merge(fgoals, by = c("lon","lat")) %>%
  merge(gfdlesm4, by = c("lon","lat")) %>%
  merge(inm, by = c("lon","lat")) %>%
  merge(inmcm5, by = c("lon","lat")) %>%
  merge(mpi, by = c("lon","lat")) %>%
  merge(mri, by = c("lon","lat")) %>%
  merge(noresm, by = c("lon","lat"))

all$cmip6.wdc <- with(all, ifelse(cat.AI.cmip6 == cat.AI.wdc, "same","diff"))
all$cmip6.era5 <- with(all, ifelse(cat.AI.cmip6 == cat.AI.era5, "same","diff"))
all$wdc.era5 <- with(all, ifelse(cat.AI.wdc == cat.AI.era5, "same","diff"))


all$casesm2.wdc <- with(all, ifelse(cat.AI.casesm2 == cat.AI.wdc, "same","diff"))
all$cesm.wdc <- with(all, ifelse(cat.AI.cesm == cat.AI.wdc, "same","diff"))
all$cmcc.wdc <- with(all, ifelse(cat.AI.cmcc == cat.AI.wdc, "same","diff"))
all$cmccesm2.wdc <- with(all, ifelse(cat.AI.cmccesm2 == cat.AI.wdc, "same","diff"))
all$cnrm.wdc <- with(all, ifelse(cat.AI.cnrm == cat.AI.wdc, "same","diff"))
all$ecearth3.wdc <- with(all, ifelse(cat.AI.ecearth3 == cat.AI.wdc, "same","diff"))
all$fgoals.wdc <- with(all, ifelse(cat.AI.fgoals == cat.AI.wdc, "same","diff"))
all$gfdlesm4.wdc <- with(all, ifelse(cat.AI.gfdlesm4 == cat.AI.wdc, "same","diff"))
all$inm.wdc <- with(all, ifelse(cat.AI.inm == cat.AI.wdc, "same","diff"))
all$inmcm5.wdc <- with(all, ifelse(cat.AI.inmcm5 == cat.AI.wdc, "same","diff"))
all$mpi.wdc <- with(all, ifelse(cat.AI.mpi == cat.AI.wdc, "same","diff"))
all$mri.wdc <- with(all, ifelse(cat.AI.mri == cat.AI.wdc, "same","diff"))
all$noresm.wdc <- with(all, ifelse(cat.AI.noresm == cat.AI.wdc, "same","diff"))

df1 <- table(all$cmip6.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMIP6-WDC"))
df2 <- table(all$casesm2.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CASESM2-WDC"))
df3 <- table(all$cesm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CESM-WDC"))
df4 <- table(all$cmcc.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMCC-WDC"))
df5 <- table(all$cmccesm2.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CMCCESM2-WDC"))
df6 <- table(all$cnrm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "CNRM-WDC"))
df7 <- table(all$ecearth3.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "ECEARTH3-WDC"))
df8 <- table(all$fgoals.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "FGOALS-WDC"))
df9 <- table(all$gfdlesm4.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "GFDLESM4-WDC"))
df10 <- table(all$inm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "INM-WDC"))
df11 <- table(all$inmcm5.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "INMCM5-WDC"))
df12 <- table(all$mpi.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "MPI-WDC"))
df13 <- table(all$mri.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "MRI-WDC"))
df14 <- table(all$noresm.wdc) %>% as.data.frame() %>% setNames(c("diffornot", "NORESM-WDC"))

df0 <- merge(df1,df2, by = "diffornot") %>% merge(df3, by = "diffornot") %>% 
  merge(df4, by = "diffornot") %>% merge(df5, by = "diffornot") %>% merge(df6, by = "diffornot") %>% 
  merge(df7, by = "diffornot") %>% merge(df8, by = "diffornot") %>% merge(df9, by = "diffornot") %>%
  merge(df10, by = "diffornot") %>% merge(df11, by = "diffornot") %>% merge(df12, by = "diffornot") %>%
  merge(df13, by = "diffornot") %>% merge(df14, by = "diffornot")

prop.diff.wdc <- df0[1,]/(df0[1,] + df0[2,])*100      
prop.diff.wdc[1,1] <- "Percent diff"
                                                      
df.prop.wdc <- rbind(df0, prop.diff.wdc)

kable(df.prop.wdc, caption = "Percent of gridcells in which the aridity category differs from Worldclim", digits = 1) %>%
  kable_styling(bootstrap_options = "bordered")

Percent of gridcells in which the aridity category differs from Worldclim

diffornot	CMIP6-WDC	CASESM2-WDC	CESM-WDC	CMCC-WDC	CMCCESM2-WDC	CNRM-WDC	ECEARTH3-WDC	FGOALS-WDC	GFDLESM4-WDC	INM-WDC	INMCM5-WDC	MPI-WDC	MRI-WDC	NORESM-WDC
diff	3273.0	4275.0	3557.0	4006.0	6503	4594.0	3691	5083.0	3531.0	4323.0	4318.0	4276.0	3349.0	3843.0
same	18418.0	17407.0	18125.0	17676.0	15179	17097.0	18000	16608.0	18160.0	17359.0	17364.0	17415.0	18333.0	17839.0
Percent diff	15.1	19.7	16.4	18.5	30	21.2	17	23.4	16.3	19.9	19.9	19.7	15.4	17.7

all$casesm2.era5 <- with(all, ifelse(cat.AI.casesm2 == cat.AI.era5, "same","diff"))
all$cesm.era5 <- with(all, ifelse(cat.AI.cesm == cat.AI.era5, "same","diff"))
all$cmcc.era5 <- with(all, ifelse(cat.AI.cmcc == cat.AI.era5, "same","diff"))
all$cmccesm2.era5 <- with(all, ifelse(cat.AI.cmccesm2 == cat.AI.era5, "same","diff"))
all$cnrm.era5 <- with(all, ifelse(cat.AI.cnrm == cat.AI.era5, "same","diff"))
all$ecearth3.era5 <- with(all, ifelse(cat.AI.ecearth3 == cat.AI.era5, "same","diff"))
all$fgoals.era5 <- with(all, ifelse(cat.AI.fgoals == cat.AI.era5, "same","diff"))
all$gfdlesm4.era5 <- with(all, ifelse(cat.AI.gfdlesm4 == cat.AI.era5, "same","diff"))
all$inm.era5 <- with(all, ifelse(cat.AI.inm == cat.AI.era5, "same","diff"))
all$inmcm5.era5 <- with(all, ifelse(cat.AI.inmcm5 == cat.AI.era5, "same","diff"))
all$mpi.era5 <- with(all, ifelse(cat.AI.mpi == cat.AI.era5, "same","diff"))
all$mri.era5 <- with(all, ifelse(cat.AI.mri == cat.AI.era5, "same","diff"))
all$noresm.era5 <- with(all, ifelse(cat.AI.noresm == cat.AI.era5, "same","diff"))

df1 <- table(all$cmip6.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMIP6-ERA5"))
df2 <- table(all$casesm2.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CASESM2-ERA5"))
df3 <- table(all$cesm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CESM-ERA5"))
df4 <- table(all$cmcc.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMCC-ERA5"))
df5 <- table(all$cmccesm2.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CMCCESM2-ERA5"))
df6 <- table(all$cnrm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "CNRM-ERA5"))
df7 <- table(all$ecearth3.era5) %>% as.data.frame() %>% setNames(c("diffornot", "ECEARTH3-ERA5"))
df8 <- table(all$fgoals.era5) %>% as.data.frame() %>% setNames(c("diffornot", "FGOALS-ERA5"))
df9 <- table(all$gfdlesm4.era5) %>% as.data.frame() %>% setNames(c("diffornot", "GFDLESM4-ERA5"))
df10 <- table(all$inm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "INM-ERA5"))
df11 <- table(all$inmcm5.era5) %>% as.data.frame() %>% setNames(c("diffornot", "INMCM5-ERA5"))
df12 <- table(all$mpi.era5) %>% as.data.frame() %>% setNames(c("diffornot", "MPI-ERA5"))
df13 <- table(all$mri.era5) %>% as.data.frame() %>% setNames(c("diffornot", "MRI-ERA5"))
df14 <- table(all$noresm.era5) %>% as.data.frame() %>% setNames(c("diffornot", "NORESM-ERA5"))

df0 <- merge(df1,df2, by = "diffornot") %>% merge(df3, by = "diffornot") %>% 
  merge(df4, by = "diffornot") %>% merge(df5, by = "diffornot") %>% merge(df6, by = "diffornot") %>% 
  merge(df7, by = "diffornot") %>% merge(df8, by = "diffornot") %>% merge(df9, by = "diffornot") %>%
  merge(df10, by = "diffornot") %>% merge(df11, by = "diffornot") %>% merge(df12, by = "diffornot") %>%
  merge(df13, by = "diffornot") %>% merge(df14, by = "diffornot")

prop.diff.era5 <- df0[1,]/(df0[1,] + df0[2,])*100      
prop.diff.era5[1,1] <- "Percent diff"
                                                      
df.prop.era5 <- rbind(df0, prop.diff.era5)   

kable(df.prop.era5, caption = "Percent of gridcells in which the aridity category differs from Worldclim", digits = 1) %>%
  kable_styling(bootstrap_options = "bordered")

Percent of gridcells in which the aridity category differs from Worldclim

diffornot	CMIP6-ERA5	CASESM2-ERA5	CESM-ERA5	CMCC-ERA5	CMCCESM2-ERA5	CNRM-ERA5	ECEARTH3-ERA5	FGOALS-ERA5	GFDLESM4-ERA5	INM-ERA5	INMCM5-ERA5	MPI-ERA5	MRI-ERA5	NORESM-ERA5
diff	3197.0	3934.0	3451.0	3819.0	6969.0	3859.0	3629.0	5402.0	3860.0	4189.0	4015.0	3633.0	2910.0	3683
same	18494.0	17748.0	18231.0	17863.0	14713.0	17832.0	18062.0	16289.0	17831.0	17493.0	17667.0	18058.0	18772.0	17999
Percent diff	14.7	18.1	15.9	17.6	32.1	17.8	16.7	24.9	17.8	19.3	18.5	16.7	13.4	17

df <- table(all$wdc.era5) %>% as.data.frame() %>% setNames(c("diffornot","ERA5-WDC"))
percent.diff <- df[1,] / (df[1,]+df[2,]) * 100
percent.diff[1,1] <- "Percent diff"
df.prop.wdc.era5 <- rbind(df, percent.diff)

kable(df.prop.wdc.era5, caption = "Percent of gridcells in which the aridity category differs between Worldclim and ERA5", digits = 1) %>% 
  kable_styling(bootstrap_options = "bordered")

Percent of gridcells in which the aridity category differs between Worldclim and ERA5

diffornot	ERA5-WDC
diff	2850.0
same	18841.0
Percent diff	13.1

4 Number of inhabitants by region

https://earthdata.nasa.gov/data/catalog/sedac-ciesin-sedac-gpwv4-popcount-r11-4.11

density <- raster("Population_density/gpw_v4_population_density_rev11_2020_1_deg.tif")
plot(density)

popcount <- raster("Population_density/gpw_v4_population_count_rev11_2020_1_deg.tif")
plot(popcount)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

popcount.df <- popcount %>%
  projectRaster(to = land_mask) %>% 
  as.data.frame(xy = T) %>% setNames(c("lon","lat","pop")) %>% 
  merge(land_mask.df, by = c("lon","lat")) %>% filter(lm == 1)

era5vApop <- merge(era5vA, popcount.df, by = c("lon","lat")) 

tab.pop.by.cat <- era5vApop %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(cat.AI, Continent) %>%
  summarise(pop = sum(pop)) 

write.csv(tab.pop.by.cat, "population.count.era5.csv")

knitr::kable(tab.pop.by.cat)

cat.AI	Continent	pop
Arid	AFRICA	9.407395e+07
Arid	ASIA	2.608739e+08
Arid	CENTRAL-AMERICA	1.102838e+07
Arid	EUROPE	3.523270e+05
Arid	EUROPE-AFRICA	3.354203e+07
Arid	NORTH-AMERICA	8.275538e+05
Arid	OCEANIA	2.483547e+05
Arid	SOUTH-AMERICA	2.714742e+06
Cold	ASIA	9.005723e+06
Cold	EUROPE	6.822776e+06
Cold	NORTH-AMERICA	8.555420e+05
Cold	POLAR	NA
Cold	SOUTH-AMERICA	6.222267e+04
Dry subhumid	AFRICA	1.131375e+08
Dry subhumid	ASIA	3.724373e+08
Dry subhumid	CENTRAL-AMERICA	2.847546e+07
Dry subhumid	EUROPE	5.195230e+06
Dry subhumid	EUROPE-AFRICA	5.933504e+07
Dry subhumid	NORTH-AMERICA	9.392606e+06
Dry subhumid	OCEANIA	7.027896e+05
Dry subhumid	SOUTH-AMERICA	1.297947e+07
Humid	AFRICA	6.854010e+08
Humid	ASIA	3.055659e+09
Humid	CENTRAL-AMERICA	1.462475e+08
Humid	EUROPE	4.999673e+08
Humid	EUROPE-AFRICA	1.452014e+08
Humid	NORTH-AMERICA	2.915502e+08
Humid	OCEANIA	1.529713e+07
Humid	POLAR	4.329376e+02
Humid	SOUTH-AMERICA	3.724919e+08
Hyperarid	AFRICA	5.328568e+07
Hyperarid	ASIA	1.970440e+07
Hyperarid	EUROPE-AFRICA	4.146658e+07
Hyperarid	SOUTH-AMERICA	5.088920e+05
Semi-arid	AFRICA	2.154724e+08
Semi-arid	ASIA	5.713649e+08
Semi-arid	CENTRAL-AMERICA	3.973955e+07
Semi-arid	EUROPE	1.841667e+06
Semi-arid	EUROPE-AFRICA	8.619193e+07
Semi-arid	NORTH-AMERICA	2.806716e+07
Semi-arid	OCEANIA	3.403064e+06
Semi-arid	SOUTH-AMERICA	2.872740e+07
NA	AFRICA	3.321493e+04
NA	ASIA	1.762974e+07
NA	CENTRAL-AMERICA	7.132541e+05
NA	EUROPE	1.481967e+05
NA	EUROPE-AFRICA	1.876606e+06
NA	NORTH-AMERICA	1.403549e+03
NA	OCEANIA	3.896168e+04
NA	POLAR	NA
NA	SOUTH-AMERICA	2.214226e+05

5 Comparison of variables datasets

vars <- c("lon","lat","Continent","Name","Acronym","tas","pr","sfcWind","Rn","ea","ET0","AI","cat.AI")
merge.vars <- c("lon","lat","Continent","Name","Acronym")
calib.wide <- select(era5vA, vars) %>%
               merge(select(wdcvA, vars), 
                     by = merge.vars) %>%
          merge(select(cmip6vA, vars), 
                by = merge.vars) %>%
          setNames(c(merge.vars,
                     "tas.era5","pr.era5", "sfcWind.era5","Rn.era5","ea.era5","ET0.era5","AI.era5","cat.AI.era5",
                     "tas.wdc","pr.wdc", "sfcWind.wdc","Rn.wdc","ea.wdc","ET0.wdc","AI.wdc","cat.AI.wdc",
                     "tas.cmip6","pr.cmip6", "sfcWind.cmip6","Rn.cmip6","ea.cmip6","ET0.cmip6","AI.cmip6","cat.AI.cmip6"))

5.1 Scatterplots with Antarctica

5.1.1 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

5.1.2 ERA5 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = "#A50026")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = "#A50026")+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

knitr::kable(data.frame(precipitation = cor(calib.wide$pr.wdc, calib.wide$pr.era5, use = "pairwise.complete.obs"),
                        temperature = cor(calib.wide$tas.wdc, calib.wide$tas.era5, use = "pairwise.complete.obs"),
                        wind = cor(calib.wide$sfcWind.wdc, calib.wide$sfcWind.era5, use = "pairwise.complete.obs"),
                        Rn = cor(calib.wide$Rn.wdc, calib.wide$Rn.era5, use = "pairwise.complete.obs"),
                        ea = cor(calib.wide$ea.wdc, calib.wide$ea.era5, use = "pairwise.complete.obs"),
                        ET0 = cor(calib.wide$ET0.wdc, calib.wide$ET0.era5, use = "pairwise.complete.obs")), 
             digits = 2, caption = "Correlation coefficients between ERA5 and Worldclim") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered")

Correlation coefficients between ERA5 and Worldclim

precipitation	temperature	wind	Rn	ea	ET0
0.89	0.99	0.68	0.67	1	0.98

ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.1.3 CMIP6 and WDC

ggarrange(

ggplot(calib.wide,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = "#FDAE61")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = "#FDAE61")+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.1.4 CMIP6 and ERA5

ggarrange(
  
ggplot(calib.wide,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = "#4575B4")+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.wide,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = "#4575B4")+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.2 Scatterplots without Antarctica

calib.noant <- calib.wide %>% filter(lat > -55)

5.2.1 Boxplots

library(ggtext)

ggarrange(
  
ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = pr.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = pr.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = pr.cmip6, fill = "CMIP6"))+
  labs(x = "Mean annual precipitation, mm", y = "", fill = "")+
  xlim(0,6000)+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = tas.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = tas.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = tas.cmip6, fill = "CMIP6"))+
  labs(x = "Surface temperature, °C", y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = sfcWind.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = sfcWind.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = sfcWind.cmip6, fill = "CMIP6"))+
  labs(x = expression("Surface wind speed, " * m~s^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = Rn.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = Rn.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = Rn.cmip6, fill = "CMIP6"))+
  labs(x = expression("Rn - G, " * MJ~m^{-2}~day^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = ea.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = ea.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = ea.cmip6, fill = "CMIP6"))+
  labs(x = "ea, kPa", y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ggplot(calib.noant)+
  geom_violin(aes(y = "Worldclim", x = ET0.wdc, fill = "Worldclim"))+
  geom_violin(aes(y = "ERA5", x = ET0.era5, fill = "ERA5"))+
  geom_violin(aes(y = "CMIP6", x = ET0.cmip6, fill = "CMIP6"))+
  labs(x = expression("ET0, " * mm~day^{-1}), y = "", fill = "")+
  scale_fill_manual(values = c("#A50026", "#FDAE61", "#DAE2F7"))+
  theme_minimal(base_size = 9),

ncol = 3, nrow = 2, common.legend = T, legend = "bottom"
)

library(SpatialPack)


coords <- calib.noant[,c("lon","lat")]

df.era5.wdc <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.wdc, pr.era5, use = "pairwise.complete.obs"), cor.test(pr.wdc, pr.era5, use = "pairwise.complete.obs")[[3]], cor(pr.wdc, pr.era5, use = "pairwise.complete.obs")^2, modified.ttest(pr.wdc, pr.era5, coords)$corr, modified.ttest(pr.wdc, pr.era5, coords)$p.value),
                        temperature = c(cor(tas.wdc, tas.era5, use = "pairwise.complete.obs"), cor.test(tas.wdc, tas.era5, use = "pairwise.complete.obs")[[3]], cor(tas.wdc, tas.era5, use = "pairwise.complete.obs")^2, modified.ttest(tas.wdc, tas.era5, coords)$corr, modified.ttest(tas.wdc, tas.era5, coords)$p.value),
                        wind = c(cor(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs"), cor.test(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs")[[3]], cor(sfcWind.wdc, sfcWind.era5, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.wdc, sfcWind.era5, coords)$corr, modified.ttest(sfcWind.wdc, sfcWind.era5, coords)$p.value),
                        Rn = c(cor(Rn.wdc, Rn.era5, use = "pairwise.complete.obs"), cor.test(Rn.wdc, Rn.era5, use = "pairwise.complete.obs")[[3]], cor(Rn.wdc, Rn.era5, use = "pairwise.complete.obs")^2, modified.ttest(Rn.wdc, Rn.era5, coords)$corr, modified.ttest(Rn.wdc, Rn.era5, coords)$p.value),
                        ea = c(cor(ea.wdc, ea.era5, use = "pairwise.complete.obs"), cor.test(ea.wdc, ea.era5, use = "pairwise.complete.obs")[[3]], cor(ea.wdc, ea.era5, use = "pairwise.complete.obs")^2, modified.ttest(ea.wdc, ea.era5, coords)$corr, modified.ttest(ea.wdc, ea.era5, coords)$p.value),
                        ET0 = c(cor(ET0.wdc, ET0.era5, use = "pairwise.complete.obs"), cor.test(ET0.wdc, ET0.era5, use = "pairwise.complete.obs")[[3]], cor(ET0.wdc, ET0.era5, use = "pairwise.complete.obs")^2, modified.ttest(ET0.wdc, ET0.era5, coords)$corr, modified.ttest(ET0.wdc, ET0.era5, coords)$p.value)))

df.era5.cmip6 <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.cmip6, pr.era5, use = "pairwise.complete.obs"), cor.test(pr.cmip6, pr.era5, use = "pairwise.complete.obs")[[3]], cor(pr.cmip6, pr.era5, use = "pairwise.complete.obs")^2, modified.ttest(pr.cmip6, pr.era5, coords)$corr, modified.ttest(pr.cmip6, pr.era5, coords)$p.value),
                        temperature = c(cor(tas.cmip6, tas.era5, use = "pairwise.complete.obs"), cor.test(tas.cmip6, tas.era5, use = "pairwise.complete.obs")[[3]], cor(tas.cmip6, tas.era5, use = "pairwise.complete.obs")^2, modified.ttest(tas.cmip6, tas.era5, coords)$corr, modified.ttest(tas.cmip6, tas.era5, coords)$p.value),
                        wind = c(cor(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs"), cor.test(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs")[[3]], cor(sfcWind.cmip6, sfcWind.era5, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.cmip6, sfcWind.era5, coords)$corr, modified.ttest(sfcWind.cmip6, sfcWind.era5, coords)$p.value),
                        Rn = c(cor(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs"), cor.test(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs")[[3]], cor(Rn.cmip6, Rn.era5, use = "pairwise.complete.obs")^2, modified.ttest(Rn.cmip6, Rn.era5, coords)$corr, modified.ttest(Rn.cmip6, Rn.era5, coords)$p.value),
                        ea = c(cor(ea.cmip6, ea.era5, use = "pairwise.complete.obs"), cor.test(ea.cmip6, ea.era5, use = "pairwise.complete.obs")[[3]], cor(ea.cmip6, ea.era5, use = "pairwise.complete.obs")^2, modified.ttest(ea.cmip6, ea.era5, coords)$corr, modified.ttest(ea.cmip6, ea.era5, coords)$p.value),
                        ET0 = c(cor(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs"), cor.test(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs")[[3]], cor(ET0.cmip6, ET0.era5, use = "pairwise.complete.obs")^2, modified.ttest(ET0.cmip6, ET0.era5, coords)$corr, modified.ttest(ET0.cmip6, ET0.era5, coords)$p.value)))

df.wdc.cmip6 <- with(calib.noant, 
                    data.frame(precipitation = c(cor(pr.cmip6, pr.wdc, use = "pairwise.complete.obs"), cor.test(pr.cmip6, pr.wdc, use = "pairwise.complete.obs")[[3]], cor(pr.cmip6, pr.wdc, use = "pairwise.complete.obs")^2, modified.ttest(pr.cmip6, pr.wdc, coords)$corr, modified.ttest(pr.cmip6, pr.wdc, coords)$p.value),
                        temperature = c(cor(tas.cmip6, tas.wdc, use = "pairwise.complete.obs"), cor.test(tas.cmip6, tas.wdc, use = "pairwise.complete.obs")[[3]], cor(tas.cmip6, tas.wdc, use = "pairwise.complete.obs")^2, modified.ttest(tas.cmip6, tas.wdc, coords)$corr, modified.ttest(tas.cmip6, tas.wdc, coords)$p.value),
                        wind = c(cor(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs"), cor.test(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs")[[3]], cor(sfcWind.cmip6, sfcWind.wdc, use = "pairwise.complete.obs")^2, modified.ttest(sfcWind.cmip6, sfcWind.wdc, coords)$corr, modified.ttest(sfcWind.cmip6, sfcWind.wdc, coords)$p.value),
                        Rn = c(cor(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs"), cor.test(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs")[[3]], cor(Rn.cmip6, Rn.wdc, use = "pairwise.complete.obs")^2, modified.ttest(Rn.cmip6, Rn.wdc, coords)$corr, modified.ttest(Rn.cmip6, Rn.wdc, coords)$p.value),
                        ea = c(cor(ea.cmip6, ea.wdc, use = "pairwise.complete.obs"), cor.test(ea.cmip6, ea.wdc, use = "pairwise.complete.obs")[[3]], cor(ea.cmip6, ea.wdc, use = "pairwise.complete.obs")^2, modified.ttest(ea.cmip6, ea.wdc, coords)$corr, modified.ttest(ea.cmip6, ea.wdc, coords)$p.value),
                        ET0 = c(cor(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs"), cor.test(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs")[[3]], cor(ET0.cmip6, ET0.wdc, use = "pairwise.complete.obs")^2, modified.ttest(ET0.cmip6, ET0.wdc, coords)$corr, modified.ttest(ET0.cmip6, ET0.wdc, coords)$p.value)))

df.combined <- rbind(
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.era5.wdc),
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.era5.cmip6),
  cbind(Metric = c("r","p-value","r2","Dutilleul","p-value"), df.wdc.cmip6)
)

knitr::kable(df.combined, 
             digits = 2, caption = "Pearson correlation coefficients and r2 between ERA5, Worldclim and CMIP6") %>%
  kableExtra::kable_styling(bootstrap_options = "bordered") %>%
  kableExtra::group_rows(group_label = "ERA5 vs WDC", start_row = 1, end_row = 5) %>%
  kableExtra::group_rows("ERA5 vs CMIP6", 6,10) %>%
  kableExtra::group_rows("WDC vs CMIP6", 11,15)

Pearson correlation coefficients and r2 between ERA5, Worldclim and CMIP6

Metric	precipitation	temperature	wind	Rn	ea	ET0
ERA5 vs WDC
r	0.88	1.00	0.77	0.89	1.00	0.97
p-value	0.00	0.00	0.00	0.00	0.00	0.00
r2	0.77	1.00	0.59	0.79	0.99	0.94
Dutilleul	0.88	1.00	0.77	0.89	1.00	0.97
p-value	0.00	0.00	0.00	0.00	0.00	0.00
ERA5 vs CMIP6
r	0.89	1.00	0.88	0.99	0.99	0.98
p-value	0.00	0.00	0.00	0.00	0.00	0.00
r2	0.79	1.00	0.77	0.97	0.98	0.96
Dutilleul	0.89	1.00	0.88	0.99	0.99	0.98
p-value	0.00	0.00	0.00	0.00	0.00	0.00
WDC vs CMIP6
r	0.89	1.00	0.78	0.91	0.99	0.96
p-value	0.00	0.00	0.00	0.00	0.00	0.00
r2	0.79	0.99	0.61	0.82	0.97	0.91
Dutilleul	0.89	1.00	0.78	0.91	0.99	0.96
p-value	0.00	0.00	0.00	0.00	0.00	0.00

5.2.2 ERA5 compared to WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.wdc, calib.noant$tas.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.wdc, calib.noant$pr.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.wdc, calib.noant$sfcWind.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.wdc, calib.noant$Rn.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC", y = "ea, kPa, ERA5", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.wdc, calib.noant$ea.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[1])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.wdc, calib.noant$ET0.era5, use = "pairwise.complete.obs"),2)), col = colorvec[1])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.3 CMIP6 and WDC

ggarrange(

ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.wdc, calib.noant$tas.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.wdc, calib.noant$pr.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.wdc, calib.noant$sfcWind.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.cmip6), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.wdc, calib.noant$Rn.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.wdc, calib.noant$ea.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[4])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.wdc, calib.noant$ET0.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[4])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.2.4 CMIP6 and ERA5

ggarrange(
  
ggplot(calib.noant,aes(x = tas.era5))+geom_point(aes( y = tas.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6")+
  lims(x = c(-30, 31), y = c(-30,31))+
   annotate(geom = "text", x =  min(c(-30, 31)) + dist(c(-30, 31))*0.15, y =  max(c(-30, 31)) - dist(c(-30, 31))*0.15, 
            label = paste("r = ", round(cor(calib.noant$tas.era5, calib.noant$tas.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = pr.era5))+geom_point(aes( y = pr.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4], colorvec[4]))+
  lims(x = c(-1, 6000), y = c(-1, 6000))+
   annotate(geom = "text", x =  min(c(-1, 6000)) + dist(c(-1, 6000))*0.15, y =  max(c(-1, 6000)) - dist(c(-1, 6000))*0.15, 
            label = paste("r = ", round(cor(calib.noant$pr.era5, calib.noant$pr.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = sfcWind.era5))+geom_point(aes( y = sfcWind.cmip6), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6")+
  lims(x = c(0,10), y = c(0,10))+
   annotate(geom = "text", x =  min(c(0,10)) + dist(c(0,10))*0.15, y =  max(c(0,10)) - dist(c(0,10))*0.15,
            label = paste("r = ", round(cor(calib.noant$sfcWind.era5, calib.noant$sfcWind.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = Rn.era5))+geom_point(aes( y = Rn.cmip6), alpha = 0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, MJ/m2/day, CMIP6", col = "Dataset")+
   lims(x = c(-8,20), y = c(-8,20))+
   annotate(geom = "text", x =  min(c(-8,20)) + dist(c(-8,20))*0.15, y =  max(c(-8,20)) - dist(c(-8,20))*0.15,
            label = paste("r = ", round(cor(calib.noant$Rn.era5, calib.noant$Rn.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ea.era5))+geom_point(aes(y = ea.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5, °C", y = "ea, kPa, CMIP6", col = "Dataset")+
   lims(x = c(0,3), y = c(0,3))+
   annotate(geom = "text",  x =  min(c(0,3)) + dist(c(0,3))*0.15, y =  max(c(0,3)) - dist(c(0,3))*0.15,
            label = paste("r = ", round(cor(calib.noant$ea.era5, calib.noant$ea.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.noant,aes(x = ET0.era5))+geom_point(aes(y = ET0.cmip6, col = "cmip6"), alpha =0.1, shape = 16, size = 0.5, col = colorvec[8])+
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
lims(x = c(-261,4200), y = c(-261,4200))+
   annotate(geom = "text",  x =  min(c(-261,4200)) + dist(c(-261,4200))*0.15, y =  max(c(-261,4200)) - dist(c(-261,4200))*0.15,
            label = paste("r = ", round(cor(calib.noant$ET0.era5, calib.noant$ET0.cmip6, use = "pairwise.complete.obs"),2)), col = colorvec[8])+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2

)

5.3 Scatterplots only Antarctica

calib.ant <- calib.wide %>% filter(lat < -60)

5.3.1 ERA5 and WDC

Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none")


ggarrange(

ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
  geom_density(aes(x=tas.era5, col = "ERA5"))+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
  scale_color_manual(values = c(colorvec[1], colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
   scale_color_manual(values = c(colorvec[1], colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),


ncol = 3, nrow = 2
)

5.3.2 CMIP6 and WDC

ggarrange(

ggplot(calib.ant,aes(x = tas.wdc))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.wdc))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.wdc))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
  scale_color_manual(values = c(colorvec[4]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.wdc))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.wdc))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.wdc))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
   scale_color_manual(values = c(colorvec[4]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2
)

5.3.3 CMIP6 and ERA5

ggarrange(
  
ggplot(calib.ant,aes(x = tas.era5))+
  geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
  scale_color_manual(values = colorvec[8])+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = pr.era5))+
  geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = sfcWind.era5))+
  geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
  geom_abline(slope = 1, intercept = 0)+
  labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
  scale_color_manual(values = c(colorvec[8]))+
  theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = Rn.era5))+
   geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
   stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ea.era5))+
   geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ggplot(calib.ant,aes(x = ET0.era5))+
   geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
  stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
   geom_abline(slope = 1, intercept = 0)+
   labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
   scale_color_manual(values = c(colorvec[8]))+
   theme_minimal()+
  theme(legend.position = "none"),

ncol = 3, nrow = 2)

5.4 Aridity category difference

5.4.1 ERA5 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.era5 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.era5, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid                Arid to Humid 
##                         1359                            1 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           18                          155 
##                 Cold to Arid                 Cold to Cold 
##                           19                        10101 
##         Cold to Dry subhumid                Cold to Humid 
##                            9                          290 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           20                           23 
##         Dry subhumid to Cold Dry subhumid to Dry subhumid 
##                            5                          205 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          209                          271 
##                Humid to Arid                Humid to Cold 
##                           30                          520 
##        Humid to Dry subhumid               Humid to Humid 
##                          347                         5436 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          210                          230 
##            Hyperarid to Cold       Hyperarid to Hyperarid 
##                            2                          726 
##                   NA to Cold            Semi-arid to Arid 
##                          570                          181 
##            Semi-arid to Cold    Semi-arid to Dry subhumid 
##                            3                          173 
##           Semi-arid to Humid       Semi-arid to Semi-arid 
##                           71                         1077
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "ERA5 compared to WDC", title = "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.era5-AI.wdc))+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), 
               palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")),
               guide = guide_legend(label.theme = element_text(angle = 0))
               )+
  labs(title = "Changes in AI", fill = "ERA5 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

ggplot() + 
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < 277 & pr.era5 - pr.wdc > - 98), aes(x=lon, y = lat,  fill = pr.era5-pr.wdc))+
  geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc > 277), aes(x=lon, y = lat), fill = "#44FF44")+
   geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < - 98), aes(x=lon, y = lat), fill = "#FF4444")+
  labs(title = "Most extreme differences in precipitation", subtitle = "In red, the 10% most extreme gricells in which Worldclim is wetter than ERA5; \nin green, the 10% most extreme gridcells in which ERA5 is wetter than Worldclim", fill = "precipitation in ERA5 - precipitations in WDC,\nmm/year")+
  theme_void(base_size = 10)+
  theme(legend.position = "bottom")

5.4.2 CMIP6 and WDC

calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                         1099                            9 
##            Arid to Hyperarid            Arid to Semi-arid 
##                           83                          150 
##                 Cold to Arid                 Cold to Cold 
##                           25                        10491 
##         Cold to Dry subhumid                Cold to Humid 
##                            5                          588 
##            Cold to Semi-arid         Dry subhumid to Arid 
##                           23                           36 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                          147                          196 
##    Dry subhumid to Hyperarid    Dry subhumid to Semi-arid 
##                            1                          350 
##                Humid to Arid                Humid to Cold 
##                           77                          138 
##        Humid to Dry subhumid               Humid to Humid 
##                          427                         5112 
##           Humid to Hyperarid           Humid to Semi-arid 
##                            1                          435 
##            Hyperarid to Arid            Hyperarid to Cold 
##                           85                            2 
##       Hyperarid to Hyperarid                   NA to Cold 
##                          652                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          520                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          111                            7 
##       Semi-arid to Semi-arid 
##                          775
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Arid to Cold", "dryer",       
                                ifelse(change.cat.wdc == "Cold to Humid", "dryer",
                                ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.wdc == "Humid to Cold","dryer",
                                ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
                                ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
                                NA))))))))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to WDC", "Change of category")+
  ylim(-55,90)+
  theme_void(base_size = 15)+
  theme(legend.position = "bottom"),

ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.wdc))+
   borders()+
  binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")), 
                                                                                                              guide = guide_legend(label.theme = element_text(angle = 0)))+
  labs(title = "Change in AI", fill = "CMIP6 compared to WDC")+
  theme_void(base_size = 15)+ylim(-55,90)+
  theme(legend.position = "bottom"),

ncol = 2)

5.4.3 CMIP6 and ERA5

calib.wide$same.cat.era5 <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.era5, "y","n"))
calib.wide$change.cat.era5 <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.era5, sep = " to "))
table(calib.wide$change.cat.era5)
## 
##                 Arid to Arid         Arid to Dry subhumid 
##                          989                            6 
##                Arid to Humid            Arid to Hyperarid 
##                            2                          230 
##            Arid to Semi-arid                 Cold to Cold 
##                          114                        10377 
##         Cold to Dry subhumid                Cold to Humid 
##                           10                          745 
##         Dry subhumid to Arid Dry subhumid to Dry subhumid 
##                           14                          169 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          219                          328 
##                Humid to Arid                Humid to Cold 
##                            5                           62 
##        Humid to Dry subhumid               Humid to Humid 
##                          382                         5456 
##           Humid to Semi-arid            Hyperarid to Arid 
##                          285                           14 
##       Hyperarid to Hyperarid                     NA to NA 
##                          725                          570 
##            Semi-arid to Arid    Semi-arid to Dry subhumid 
##                          511                          146 
##           Semi-arid to Humid       Semi-arid to Hyperarid 
##                          121                            3 
##       Semi-arid to Semi-arid 
##                          778
calib.wide$dryer.era5 <- with(calib.wide, ifelse(change.cat.era5 == "Arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Arid to Semi-arid", "dryer",
                                ifelse(change.cat.era5 == " Arid to Dry subhumid", "dryer", 
                                ifelse(change.cat.era5 == "Arid to Hyperarid", "wetter",
                                ifelse(change.cat.era5 == "Cold to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Cold", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Dry subhumid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Dry subhumid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Semi-arid", "wetter",
                                ifelse(change.cat.era5 == "Humid to Cold","dryer",
                                ifelse(change.cat.era5 == "Humid to Dry subhumid","wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Dry subhumid", "dryer",
                                ifelse(change.cat.era5 == "Semi-arid to Arid", "wetter",
                                ifelse(change.cat.era5 == "Semi-arid to Humid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Arid", "dryer",
                                ifelse(change.cat.era5 == "Hyperarid to Semi-arid", "dryer",
                                NA)))))))))))))))))))

ggarrange(
ggplot(filter(calib.wide, same.cat.era5 == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.era5))+
  scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
  borders()+
  labs(fill = "CMIP6 compared to ERA5", "Change of aridity category")+
  ylim(-55,90)+
  theme_void()+
  theme(legend.position = "bottom"),


ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat,  fill = AI.cmip6-AI.era5))+
   borders()+
  binned_scale(aesthetics = "fill", breaks =c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")))+
  labs(title = "Change of AI", fill = "CMIP6 compared to ERA5")+
  theme_void()+ylim(-55,90)+
  theme(legend.position = "bottom"),
ncol = 2)

6 Aridity Index Hargreaves

calib$

7 NDVI

# Download Gimms NDVI
library(gimms)
library(gganimate)

land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land

#ndvi.ref <- downloadGimms(x= 1981,y= 2000,dsn = "NDVI")

list.ndvi <- list.files(path = "NDVI", pattern = "ndvi.*\\.nc4$", full.names = T)[1:2]

ndvi.rast <- rasterizeGimms(list.ndvi, keep = 0) %>% mean(na.rm = T) %>% projectRaster(land_mask) # rasterize Gimms does the quality control and removes NAs
ndvi.df <- ndvi.rast %>% as.data.frame(xy = T) %>% setNames(c("lon","lat","ndvi")) %>% merge(land_mask.df, by = c("lon","lat")) %>% filter(lm != 0)

write.table(ndvi.df, "NDVI/ndvi.df.txt")

ndvi.df <- read.table("NDVI/ndvi.df.txt") 

ggplot(subset(ndvi.df, lm == 1))+geom_raster(aes(x=lon, y = lat, fill = ndvi))+
  scale_fill_continuous_sequential(palette = "Greens 3", rev = T)+
  ylim(-55,90)+labs(fill = "NDVI")+
  theme_void()+theme(legend.position = "bottom")

calibn <- merge(calib, ndvi.df, by = c("lon","lat"))

ggplot(subset(calibn, !is.na(cat.AI)))+geom_boxplot(aes(x=cat.AI, y = ndvi, col = cat.AI))+facet_grid(rows = vars(source))+
  scale_color_manual(values = col.cat)+
  labs(x = "", y = "NDVI")+
  theme_minimal()+theme(legend.position = "none")

aov(ndvi~cat.AI, data = group_by(calibn, source)) %>% summary()
##                Df Sum Sq Mean Sq F value Pr(>F)    
## cat.AI          5   1026  205.17    8708 <2e-16 ***
## Residuals   43200   1018    0.02                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23577 observations effacées parce que manquantes

range.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(range.ndvi = paste(round(quantile(ndvi, na.rm = T)[c(2,4)],2), collapse = " to ")) %>%
  cast(source~cat.AI)

cor.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>% 
  summarise(pearson.r = round(cor(ndvi, AI, use = "pairwise.complete"),2)) %>%
  cast(source~cat.AI)

knitr::kable(range.ndvi, caption = "NDVI range by source and aridity category")

NDVI range by source and aridity category

source	Cold	Hyperarid	Arid	Semi-arid	Dry subhumid	Humid
CMIP6	0.3 to 0.47	0.09 to 0.11	0.1 to 0.19	0.17 to 0.36	0.24 to 0.49	0.4 to 0.69
ERA5	0.27 to 0.47	0.09 to 0.11	0.11 to 0.21	0.19 to 0.37	0.24 to 0.48	0.42 to 0.68
Worldclim	0.29 to 0.47	0.09 to 0.11	0.1 to 0.19	0.18 to 0.34	0.28 to 0.44	0.45 to 0.7

knitr::kable(cor.ndvi, caption = "Correlation between NDVI and AI by source and aridity category")

Correlation between NDVI and AI by source and aridity category

source	Cold	Hyperarid	Arid	Semi-arid	Dry subhumid	Humid
CMIP6	-0.08	0.15	0.36	0.38	-0.01	0.31
ERA5	-0.16	0.08	0.47	0.30	0.08	0.38
Worldclim	0.05	0.05	0.49	0.39	0.14	0.52

knitr::knit_exit()

7.0.1 Maps CMIP6 ERA5 // WDC

In red: when CMIP6 or ERA5 are warmer than Wordlclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.era5-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: ERA5 - WDC")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - WDC")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or ERA5 are wetter than Worldclim

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or ERA5 compared to Worldclim

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.era5-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: ERA5 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.wdc))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - WDC", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

7.0.2 Maps CMIP6 WDC // ERA5

In red: when CMIP6 or WDC are warmer than ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.wdc-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: WDC - ERA5")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - ERA5")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when CMIP6 or WDC are wetter than ERA5

cowplot::plot_grid(
  
  ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.wdc-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)

In blue: when wind speed is higher for CMIP6 or WDC compared to ERA5

cowplot::plot_grid(
  
ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.wdc-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: WDC - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom"),

ggplot(calib.wide)+
  geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.era5))+
   scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - ERA5", guide = "legend")+
  theme_void()+
  theme(legend.position = "bottom")
)